import numpy as np
import Mixing_rules as MR
from Pure_VLE_TB import T_BOIL_CAL
from EOS_collection import CUBIC_EOS
from scipy.optimize import minimize

R = 8.314

#输入对应态参数
Tc = np.array([[305.32, 539]])
Pc = np.array([[4.872e6, 2.74e6]])
w = np.array([[0.09949, 0.351]])
Z = np.array([[0.279, 0.261]])
# M = np.array([[44, 44]]) #分子量，暂时用不上
# Tb = np.array([[322.02, 373.42]])

#输入要计算的状态参数

z1 = np.array([0.2654])
z2 = np.array([0.7346])
z = np.concatenate([z1, z2])
z = z.reshape(-1, 1)
x = z
y = z
T = np.array([430])
P = 10e5
LtoF = 0.5

EOS = CUBIC_EOS(Pc, Tc, Z, w)

def Fraction(LtoF, x, y, z, P, T):
    # 计算液相逸度系数
    sumxy = 1
    LtoF_new = LtoF
    a = EOS.PR_a(T)
    b = EOS.PR_b()

    while abs(sumxy) > 1e-8:
        LtoF = LtoF_new
        aL = MR.a_mix(a, x)
        bL = MR.b_mix(b, x)
        VL = 2 * bL
        VLnew = EOS.PR_V(P, VL, T, aL, bL)
        while np.sum ( np.abs(VLnew - VL) ) / np.sum ( VL ) > 1e-5:
            VL = VLnew
            VLnew = EOS.PR_V(P, VL, T, aL, bL)
        ZML = P * VL / R / T
        xa = MR.xya_mix(a, x)
        phiL = EOS.PR_PHI_MIX(P, VL, T, aL, bL, xa, b, ZML)
        # print(f"phiL={phiL}")

        # 计算气相逸度系数
        aV = MR.a_mix(a, y)
        bV = MR.b_mix(b, y)
        VV = R * T / P
        VVnew = EOS.PR_V(P, VV, T, aV, bV)
        while np.sum ( np.abs(VVnew - VV) ) / np.sum ( VV ) > 1e-5:
            VV = VVnew
            VVnew = EOS.PR_V(P, VV, T, aV, bV)
        ZMV = P * VV / R / T
        ya = MR.xya_mix(a, y)
        phiV = EOS.PR_PHI_MIX(P, VV, T, aV, bV, ya, b, ZMV)
        # print(f"phiV={phiV}")

        k = phiL / phiV
        x = z / (k.reshape(-1, 1) + LtoF * (1 - k.reshape(-1, 1)))
        y = k.reshape(-1, 1) * x

        sumxy = np.sum((k - 1) * z.flatten() / (k + LtoF * (1 - k)))
        term1 = sumxy
        term2 = np.sum(z.flatten() * (k - 1) ** 2 / (k + (1 - k) * LtoF) ** 2)
        LtoF_new = LtoF - term1 / term2
        # print(f"term1 = {sumxy}")
        # print(f"term2 = {term2}")

    print(f"液相组成 = {x}")
    print(f"气相组成 = {y}")
    return LtoF_new

LtoF = Fraction(LtoF, x, y, z, P, T)
print(f"液相分率是{LtoF}")